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We examine in detail the theoretical underpinnings of previous successful applications of local 
molecular field (LMF) theory to charged systems. LMF theory generally accounts for the averaged 
effects of long-ranged components of the intermolecular interactions by using an effective or restruc- 
tured external field. The derivation starts from the exact Yvon-Born Green hierarchy and shows 
that the approximation can be very accurate when the interactions averaged over are slowly varying 
at characteristic nearest-neighbor distances. Application of LMF theory to Coulomb interactions 
alone allows for great simplifications of the governing equations. LMF theory then reduces to a sin- 
gle equation for a restructured electrostatic potential that satisfies Poisson's equation defined with a 
smoothed charge density. Because of this charge smoothing by a Gaussian of width a, this equation 
may be solved more simply than the detailed simulation geometry might suggest. Proper choice 
of the smoothing length a plays a major role in ensuring the accuracy of this approximation. We 
examine the results of a basic confinement of water between corrugated wall and justify the simple 
LMF equation used in a previous publication. We further generalize these results to confinements 
that include fixed charges in order to demonstrate the broader impact of charge smoothing by a. 
The slowly-varying part of the restructured electrostatic potential will be more symmetric than the 
local details of confinements. 



I. INTRODUCTION 

The treatment of long-ranged Coulomb interactions re- 
mains a substantial challenge in nearly all classical molec- 
ular simulations. The basic problem is that interactions 
from even very distant charges remain important, as indi- 
cated by the divergence of the integral of 1 /r from any fi- 
nite tnmcation distance to infinity. This causes problems 
in computer simulations of systems with charges, where 
contributions from distant periodic images of the simula- 
tion cell must be taken into account, typically by Ewald 
sum techniques. This also hampers the development of 
simple, intuitive local pictures and analytic theory when 
charges are present. 

Here we discuss basic features of a new approach 
for charged systems, local molecular field (LMF) the- 
ory. LMF theory was originally applied to Lennard- 
Jones (LJ) systems for both simulation and analytical 
work [ll, [2, S, i, [1], and the LMF approach has 
been used recently to analyze many charged systems, 
both uniform and nonuniform, with even greater suc- 
cess 0,ii[lS[ll|. 

LMF theory can be generally characterized as a map- 
ping that relates structural and thermodynamic proper- 
ties of a nonuniform system with long-ranged and slowly- 
varying intermolecular interactions in a given external 
potential energy field, representing, e.g., fixed solutes or 
confining walls, to those of a simpler "mimic system" 
with properly chosen short-ranged intermolecular inter- 
actions in an effective or restructured field. The strong 
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short-ranged interactions in the mimic system are chosen 
to give an accurate description of the forces between typi- 
cal nearest-neighbor molecules in the full system, e.g., the 
repulsive molecular cores in a dense LJ fluid, local hy- 
drogen bonds in water, or ion pairs in a strongly coupled 
ionic system. Because of the short range of the interac- 
tions in the mimic system, fast and efficient simulations 
scaling linearly with system size are possible. The re- 
structured field is determined self-consistently and is the 
sum of the bare external field in the original system and 
a density-weighted mean-field average of the remaining 
long-ranged and slowly-varying parts of the intermolec- 
ular interactions. The restructured field corrects major 
errors that can arise, especially in nonuniform systems, 
from applying solely a spherical truncation of the inter- 
molecular pair interactions, exemplified by shifted force 
truncations for nonuniform Coulomb systems (l^ . [l^. 
More advanced spherical truncations of 1/r exist, such 
as site-site reaction field techniques , but the difficul- 
ties in nonuniform situations remain the same. 

The LMF approach offers a very general perspective 
and the averaging can be carried out for any slowly- 
varying component of the intermolecular interactions. 
We show here that LMF theory takes an especially pow- 
erful and simple form, with strong analogies to classical 
electrostatics, when it is applied uniformly to the basic 
Coulomb 1/r interaction alone (and not, e.g., to LJ inter- 
actions as well) . We can then take advantage of symme- 
tries between charges of different signs and magnitudes, 
interacting with the same spherically symmetric 1 /r po- 
tential, and determine the effective field using only the 
total charge density, and not separate number densities 
for each charged species or interaction site as would be 
required for a general mixture. We find that the con- 
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tribution to the restructured field in LMF theory from 
the slowly-varying parts of the Coulomb interactions can 
be exactly determined from a restructured electrostatic 
potential that satisfied Poisson's equation, but with a 
Gaussian-smoothed charge distribution. A proper choice 
of the smoothing length a in the Gaussian plays a key 
role in the accuracy of LMF theory, since it also deter- 
mines the size of the local strong-coupling region within 
which strong forces from nearest-neighbor sites must be 
accurately captured by the short-ranged interactions in 
the mimic system. 

Previous discussions of LMF theory for Coulomb or LJ 
systems relied on general arguments and did not exploit 
all the important simplifications that arise by applying 
LMF theory only to 1/r with the same potential split- 
ting regardless of point charge details. Here we seek to 
justify and document LMF theory as completely as pos- 
sible, accounting for Coulomb symmetries. We show how 
this leads to a simple, broadly applicable equation for a 
restructured electrostatic potential, and then we extend 
results from [ll| to show that smooth LMF solutions can 
still arise from a general class of molecularly-corrugated 
confinement potentials. 



II. DERIVATION OF THE LMF EQUATION 

The LMF equation may be derived for mixtures and 
standard site-site molecular models. However the ba- 
sic features of the derivation are identical regardless of 
whether we consider a mixture, a site-site molecular fluid, 
or a single-component system. Therefore, we present 
the derivation for a single-component system in order to 
highlight the physical approximations in the derivation 
without the clutter of the multiple indices necessary in 
the derivations for mixtures or molecules. Furthermore, 
symmetries available to us when treating charge inter- 
actions will allow us to simplify the LMF equation for 
mixtures into a single electrostatic potential equation, as 
we will show later in section Hill 



A. Motivation 

We assume that the intermolecular interactions w{r) 
in the full system of interest are slowly varying at large 
separations. The first step in LMF theory is to divide 
w(r) into properly chosen short- and long-ranged parts: 

w{r) = U(3{r) + ui{r), (1) 

such that ui(r) contains all the the slowly- varying long- 
ranged interactions in u;(r) but remains slowly vary- 
ing at characteristic nearest-neighbor distances where 
there are strong forces between neighboring molecules. 
By construction, the short-ranged component uo(r) then 
captures those strong short-ranged forces but vanishes 
quickly at larger separations. 



LMF theory then defines a mapping from the full sys- 
tem with pair potential w{r) and an external single- 
particle potential energy function (j){v) to a mimic sys- 
tem defined by the short-ranged pair interactions uo(r) 
and an effective or restructured external potential energy 
function 0R(r): 

Full Mimic 
J w{r) \ LMF J wo(r) 1 (2) 
\ 0(r) / ^ \ 0R(r) / • 

We call the final combination a mimic system when 
uo{r) and (/)r are properly chosen because we seek a sys- 
tem that can capture many relevant structural and ther- 
modynamic properties of the full system. In addition, 
we seek this mapping to be closed, in the sense that the 
mapping should not require any information about the 
full system other than the original w{r) and 4>{v). It 
seems plausible that some effective could account for 
averaged effects of the neglected long-ranged interactions 
ui{r) in the mimic system, but we seek a more careful 
development of this intuition. 

B. Exact Starting Point for the LMF Derivation 

The derivation of the LMF equation starts with an 
exact statistical mechanical equation involving structure 
and forces, the Yvon-Born-Green (YBG) hierarchy of 
equations [isl . [l6j . In one form, the first equation in 
the hierarchy reads 

fcBTVlnp(^) (r) = -V(/)(r) - j dr' p (r'|r) Vw (|r - r'|) . 

(3) 

This equation relates the gradient of the (singlet) den- 
sity profile of the system, p^^^ (r), to the force from the 
external single particle potential, (/>(r), and the averaged 
force between particles, weighted by the conditional sin- 
glet density, p(r'|r). This function p(r'|r) is defined as 
the density at position r' given that a particle is at r, 
and may be expressed mathematically as the ratio of the 
two-particle density at r and r' and the single-particle 
density at r: 



Thus, the YBG hierarchy expresses the rapidly-varying 
singlet density in terms of the more complex and even 
more rapidly-varying two-particle density. There are two 
related difficulties in trying to extract practically useful 
information about the singlet density from ([3]). First, the 
density profile is related to more complicated rapidly- 
varying functions. Second, there is no obvious manner 
to create a self-contained loop to determine the singlet 
density profile, again due to its relation to the higher- 
order density profile. 

Typical superposition closures of the YBG hierarchy 
attempt to resolve both difficulties by approximating 



3 



p(r'|r) ~ p^^") (r'). leading to an equation where the gra- 
dient over r can be pulled out of ^ and p^^^ (r) is de- 
fined in terms of itself. However, even for a bulk fluid, 
this is a poor approximation for any moderately dense 
liquid fvA [l8| . The superposition approximation ne- 
glects any excluded volume or specific-binding interac- 
tions that alter the small |r — r'| nature of p (r'|r) rela- 
tive to p^^^ (r'). Replacing the conditional singlet density 
with the singlet density simplifies the equation but also 
effectively includes with large weight too many configu- 
rations where two particles are improbably close to each 
other, leading to great inaccuracies. 

We assume that there exists a good choice of uo{r) 
that captures those short-ranged excluded-volume and 
specific-binding contributions. We will describe such a 
choice for splitting l/r using the smoothing length a 
in section Hill Given uo{r), we write the YBG equations 
for both the full sytem and the LMF-mapped system, 
denoted by the subscript R, from ^ as: 

fcBTVlnp^i) (r; [0]) ^ -V0(r) 

- I drV(r'|r; [<^]) (|r - r'|) , 

fceTVlnpl,^' (r; [0r]) = -V0R(r) 

-J dr'pR(r'|r;[0R])Vuo(|r-r'|). (5) 

The density profiles and conditional singlet densities are 
written as functionals of 4>(r) and 0R(r) to emphasize 
the restructuring of the external potential energy in the 
LMF mapping. 

As written with an arbitrary choice of uo{r), ([5]) simply 
displays formally exact but not particularly useful equa- 
tions for two arbitrary and unconnected systems with 
different intcrmolccular interactions in different external 
fields. As in we now connect these systems and 
achieve substantial simplifications and cancellations by 
requiring that (/)r be chosen so that 

pW(r;M) = pW(r;[0R]). (6) 

It seems very plausible that such a 0r exists, especially 
when uo{r) captures the strong short-ranged intermolec- 
ular forces. We then take the difference between the 
equations in ([5]). The terms involving the singlet den- 
sities on the left-hand sides cancel and the result can be 
exactly written in a form useful for further analysis as 

- V0R(r) = -V0(r) - J rfr'pR(r'; [0r])Vmi (|r - r'|) 

(7a) 

- J dr' [p(r'|r; [0]) - pR(r'|r; [<Ar])] V^^o (|r - r'|) 

(7b) 

- j dr' [p(r'|r; [0]) ~ p(r'; [0])] V^^i (|r - r'|) . 

(7c) 



Henceforth, for simplicity, we denote p'^^^ (r) as p(r). 

Equation (7) formally holds for any choice of UQ{r) 
and ui{r), and because conditional densities still appear, 
it generally has most of the problems of the standard 
hierarchy equations noted after ([?]). Thus it may seem 
we have made little progress. But this is not the case 
for proper choice of uo(f) and ui{r) with the properties 
sketched above. This will allow considerable simplifica- 
tion of (7) and defines a mimic system that 

• has a wcll-choscn uo(f) that mimics higher order 
correlations in the full system as well as the singlet 
density profile, 

• can be treated via statistical mechanics without 
any a priori knowledge of the full system aside from 
(/)(r) and w(r), and 

• has an associated 0R(r) that depends only on sin- 
glet density profiles of the mimic system and not 
on the more complex conditional singlet densities. 

C. Approximations to Yield LMF Equation 

Physically-motivated approximations to the exact (7) 
will result in a final LMF equation below in ([S]) that has 
a simple mean-field form. The following discussion elu- 
cidates that the LMF equation is not a blind assertion 
of mean-field behavior but rather a controlled and accu- 
rate approximation, provided that we choose our mimic 
system carefully. 

Recalling from ([1]) that w{r) = UQ{r) -f ui(r), the 
LMF derivation focuses on making two connected and 
physically-reasonable approximations based on choosing 
a short-ranged u^ir) that will induce the correct nearest- 
neighbor structure, where the conditional singlet and sin- 
glet densities differ the most, and a corresponding ui(r) 
that is slowly-varying on that length scale. As we will 
see in section IHI Al the form of the Coulomb potential 
grants us the freedom to fulfill both conditions by choos- 
ing a single smoothing length a. When uo(f) and ui(r) 
arc chosen correctly, we shall have a "mimic system." 

By making approximations to (7), the equality of den- 
sities expressed in ([6|) will only hold approximately, and 
we will find an equation based on the exact (7) that will 
also no longer be exact. But all will be accurately satis- 
fied with proper choices of uo{r) and ui(r). The general 
arguments below have been briefly sketched before 0], 
but the form of (7) allows us to make a more precise 
statement of the necessary approximations. 

We note that ((7a|) depends only on the density response 
of the mimic system, and no information about the den- 
sity response of the full system or any conditional density 
is needed. Thus if terms (|7b)) and ((7c)) could be neglected, 
then ([7a|) alone would define a closed equation for the 
mimic system. The right side of (j7ap still transforms as a 
gradient and should display all the associated properties. 
In that sense, neglecting terms (j7bp and (|7c)) together 
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may well be more accurate than taking cither one of the 
following approximations individually. 

With proper choice ofuo{r) and ui{r), we may reason- 
ably neglect the final two terms in (7) using the following 
arguments: 

1. Term (j7bp probes the difference between the condi- 
tional singlet densities for the full and mimic sys- 
tems via convolution with Vito(r). The integrand 
will be quickly forced to zero at larger |r — r'| by the 
vanishing gradient of the short-ranged uq (|r — r'|). 
Since both the full and mimic systems have the 
same strong short-ranged core forces with an ap- 
propriately chosen ito(r), and these core forces 
should mainly determine the short-ranged part of 
the conditional densities, it seems plausible that 
with proper choice of uo(r), term ()7bp can be ne- 
glected. This is significantly better than any super- 
position approximation for the conditional densities 
and does not require us to numerically estimate ei- 
ther conditional density. 

2. Term ([7c)) probes the difference between the con- 
ditional singlet density and the singlet density of 
the full system. As explained previously with rela- 
tion to standard closures of the YBG equation, as- 
suming their equality can be highly problematic at 
short distances. However, we are saved by the fact 
that this difference is paired with Vui(|r— r'|). 
Since uo(r) has been chosen to encompass core 
interactions, ui{r) will be simultaneously slowly- 
varying over those nearest-neighbor distances, so 
the associated force is essentially zero for exactly 
the range of distances where the conditional sin- 
glet density and singlet density differ significantly. 
Thus we expect for many liquids that the integrand 
in term (|7cp may also be accurately approximated 
as zero. 

These two approximations allow us to truncate the YBG 
hierarchy. No explicit knowledge of the conditional sin- 
glet density in either the full or the mimic system will be 
required, but we expect the conditional singlet density 
of the mimic system to closely resemble that of the full 
system at short range. 

Notably, the second approximation will fail when there 
are long-ranged pair correlations arising from capillary 
waves at the liquid-vapor interface or the divergence of 
the correlation length near the critical point. In those 
instances, the conditional singlet density will not be ap- 
proximately the same as the singlet density well beyond 
nearest-neighbor distances. However, we do expect this 
approximation to hold in strongly-coupled charged sys- 
tems away from the critical point where there exist pair 
correlations that decay exponentially over those nearest- 
neighbor distances. 

With these two coupled approximations, only (|7a|) re- 
mains, and we can take the gradient with respect to r 
outside of the integral. Then, integrating over r yields 



the general LMF equation 

Mr) = 0(r) + I dv'p^iv'; [cj,^])u, (|r - r'|) + C, (8) 

where the integration constant C will be set by bound- 
ary conditions. The form of this equation highlights 
the connection with mean-field techniques. The func- 
tion may be viewed as an averaging of the slowly- 
varying, long-ranged portions of the pair potential over 
the single-particle density /5R(r). The hierarchy has been 
successfully truncated since is defined in terms of only 
PR(r). Provided sufficient computer time to simulate the 
short-ranged mimic system and therefore close the self- 
consistent loop between 0R(r) and /OR(r), the LMF equa- 
tion is exactly soluble. 

This equation is not a naive mean-field ansatz for (/)r 
since it arises from reasonable approximations to exact 
statistical mechanical equations. The crux of LMF the- 
ory, however, is choosing reasonable uoijr) and ui(r) such 
that the necessary approximations are accurate. The 
next section discusses this choice for 1/r. 



III. ELECTROSTATICS VIA LMF THEORY 

Here, we explain our choice of UQ{r) and iti(r) for 
charge interactions. Then, we seek to highlight the sym- 
metry and simplifications possible when applying the 
LMF equation to purely charge-charge interactions by 
developing the LMF equation for a single restructured 
electrostatic potential, rather than for separate potential 
energy functions for each atom type. 



A. Coulomb Potential Split 

In CGS units, the pair interaction between two 
charges, qa and g-y, is 

wi^r) = Ml. (9) 
er 

The factor e is simply the dielectric constant of a uni- 
form medium in which the charges are immersed. The 
details of charge magnitudes, signs, and choice of CGS or 
SI units will change from system to system, but the un- 
derlying functional form 1/r will remain constant. Thus, 
following p^ . we introduce Va{r) and to split this 

function as 

i = z;o(r)+wi(r), (10) 
r 

into short-ranged and long-ranged components. Both 
VQ{r) and vi(r) will yield either a potential energy or 
an electrostatic potential, with appropriate prefactors. 

Work of previous researchers 01 revealed that a benefi- 
cial choice of vi (r) is the functional form associated with 
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FIG. 1: Demonstration of 1/r potential split for two different 
values of the smoothing length a. When o is increased from 
1 to 4, more of the 1/r core interaction is included in voir), 
and v\(r) becomes correspondingly more slowly varying. 



the electrostatic potential due to a unit Gaussian charge 
distribution of width a, defined as 



r3/2tT' 



■ exp 



leading to a vi (r) given by 
vi{r) = / dv'pciv') 



erf (?-/ct) 



r - r' 



(11) 



(12) 



This vi{r) is slowly- varying in r-space over the smooth- 
ing length a. As discussed in 0], this choice of ui(r) also 
simultaneously isolates only the small wave-vector com- 
ponents of the Coulomb potential. The Fourier transform 
of v\(r) is 



(13) 



and the Gaussian function attenuates any large-A: compo- 
nents that would be rapidly- varying in r-space. Thus the 
wi(r) is particularly well-suited for mean-field averaging. 

The Wo (r) corresponding to our choice of v\ (r) will then 
be given as 



wo(f) 



erfc(r/cr) 



(14) 



Therefore, v^ir) = \ / r — vi[r) is proportional to the elec- 
trostatic potential due to a point charge screened by a 
surrounding, neutralizing Gaussian charge distribution 
whose width a also sets the scale for the smooth trunca- 
tion of vq (r) . This split is shown in figure [T] along with 
a demonstration of the impact of increasing the smooth- 
ing length <j on wo(r) and ui(r). The form of fi(r) is 
slowly varying over ct, and VQ{r) approximately repre- 
sents a Coulomb core potential of range a. Tuning a to 
larger values increases the range of VQ{r) and simultane- 
ously smooths the variation of Ui(r) over characteristic 
nearest- neighbor interactions. 



Since 1/r has no inherent length scale, the split be- 
tween short- and long-ranged by the smoothing length a 
may be tailored to each situation. This length a should 
not be viewed as a fitting parameter, but rather a consis- 
tency parameter. We expect that choosing a too small 
will result in poor and rapidly-varying results from full 
LMF theory as a is changed since the mean-field average 
will be over rapidly-varying interactions. However once 
a reaches a (Tmin that is related to characteristic nearest- 
neighbor distances, the assumptions used in the previous 
section to derive the LMF equation from the YBG hier- 
archy become accurate, and essentially the same results 
should be found for larger a values as well [l^ . Thus we 
expect that self-consistently applied LMF theory with an 
appropriately chosen a will yield a true mimic system. 
This requirement for reasonable core inclusion has been 
explored in previous papers 0, B @ • 

For small r < cr, the force due to VQ{r) is quite similar 
to that of a bare point charge. By increasing a, we in- 
crease the effective range of these essentially unscreened 
Coulomb interactions. Thus, v^ir) can be thought of as 
a Gaussian-truncated (GT) or "short" Coulomb inter- 
action. We refer to a model with generically-truncated 
Coulomb interactions as a short model and with VQ{r) 
specifically used as a GT model. We showed in [lH that 
short GT water, where all the point charges in SPC/E 
water are replaced by VQ{r) with no account taken of 
the effective field, can give a surprisingly accurate de- 
scription of all site-site correlations for bulk water. But 
this same approach can fail spectacularly for electrostatic 
correlations in nonuniform systems unless corrected by 
the inclusion of the l ong -ranged forces via the LMF- 
restructured potential [l]| . 

The functional form of vn{r) suggests both Ewald 
sums [1^ and Wolf sums [2l|. Ewald sums, however, 
seek to explicitly evaluate the forces from vi{r) at each 
timestep using an exact sum over the fluctuating peri- 
odic images of the simulation cell, and thus have a very 
different philosophical approach. The analogy with Wolf 
sums is much more appropriate in that each results in 
a spherical truncation vo{r) applied at each time step. 
However, the historical development of the Wolf summa- 
tion technique relies on the observation that charges in 
ionic lattices are quickly neutralized by slightly shifted 
"mirror" lattices of neutralizing charges. In fluids, this 
is altered to assuming that charges within a sharp cutoff 
radius are exactly neutralized by a shell of charge at the 
cutoff radius, and the use of erfc(r/<T) to damp the forces 
from 1/r is applied more as a numerical afterthought. 
In the context of LMF theory, we view such damping as 
much more crucial to the success of Wolf sums in repro- 
ducing the structure of uniform fluids than the imposition 
of strict neutrality within a sharp cutoff radius. 

The GT vo{r) also bears great similarity to site-site 
reaction field (RF) approaches [13, [13, [H, [11], when 
the fluid is assumed to be conducting. The usual RF 
method generates a truncated i'o,rf('") from the electro- 
static potential of a unit point charge surrounded by a 
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uniform neutralizing spherical charge distribution. We 
discuss in the Appendix more details of the relation be- 
tween the GT vo{r), and the usual RF vo,rf('') as well as 
a smoother charged-clouds RF truncation 'yo,cc('') Hi- 
The choice of spheres as the neutralizing charge distri- 
bution in these RF methods leads to a potential that is 
exactly zero beyond a certain cutoff distance, but this in- 
duces a discontinuity in a higher order derivative at the 
sharp cutoff. The LMF choice of Gaussians leads to an 
infinitely smooth and soft truncation both in r-space and 
/c-space. Moreover a GT model logically separates the 
choice of the smoothing length a from the possible use for 
computational convenience of a sharp cutoff at a larger 
distance along with various local smoothing schemes that 
could be employed near the cutoff. 

Despite the general similarities, the theoretical con- 
structs of both reaction field techniques and Wolf sums 
allow for no obvious correction in nonuniform systems, 
such as the slab confinement to which we applied LMF 
theory in [llj. As such, LMF theory is an important 
advance in the use of spherical truncations for l/r. 

B. LMF Electrostatic Potential 

Now we will derive a further simplified electrostatic 
form of the LMF equation. We need only consider the 
LMF mixture equation here because, with appropriate 
further approximations, the LMF equation for standard 
molecular models collapses onto the mixture LMF equa- 
tion as will be shown in a paper dealing more broadly 
with site-site molecular models [l^. Using the YBG 
equation for mixtures, exactly the same approximations 
as in section [n] lead to the following LMF equation for a 
species a in a mixture of different species indexed by 7 
and including a, as originally given in Q: 

(l)R,a = <t)air) + dr' pR^^{r') ■ ui^at{\r - r'\). (15) 

If the pair interactions being treated were of general form 
this would be the furthest simplification possible, with an 
LMF equation to be solved for each different species. 

The observation in Q that a single a for all charged 
interactions proves most useful leads straightforwardly to 
a much more compact form. Using the LMF approach 
only for charges, each ui,Q'y(r) may be written as 

ggg^ erf (r/o-) 
wi,Q7(f) = — ^ = —^vi{r). (16) 

Thus (|15p may be exactly written as 

'^R,a(i-) = 0a(i-) + y y" dr' ^^g^pR,^(r')j -widr - r'l). 

(17) 

Using the natural definition of charge density as 

p'(r)=^g^p^(r), (18) 
7 



we may reexprcss the previous equation for (/>R,,a as 

0R,a(r) = 0.(r) + ^ I drV^(r')•^.l(|r-r'|). (19) 

By noting as done in Q that may be due to both 
electrostatic and nonelectrostatic components, we sepa- 
rate (f>a as 

0a(r) = (/)nc,a(r)-f gaV(r) (20) 

where (f)nc,a encompases nonelectrostatic confinements 
such as the smooth LJ walls used in [llj and V repre- 
sents the general external electrostatic potential. 

This observation naturally leads to the following 
rewriting of each LMF equation for a species a as 

0R,a(r) = 0nc,Q(r) + gaVR(r), (21) 

with a single restructured electrostatic potential Vr, de- 
fined as 

VR(r) = V(r) + i 1 drV^(r') • «i(|r - r'|). (22) 

Thus all self-consistent requirements on the diverse set 
of ^R,a can be written in terms of a single restructured 
electrostatic potential Vr defined by one LMF equation 
involving the equilibrium mobile charge density p"^ . 

Given the convolution definition of vi{r) as the po- 
tential due to a Gaussian charge density in we 
may schematically represent the integral term in (|22p as 
p'^*PG*l/r where * indicates a convolution. Thus p2|) 
can be exactly rewritten as 

VR(r) = V(r) + -J rfrV|"(r') • (23) 

where p'" represents the smoothing of the charge density 
by pq in (jlip via the convolution 

p^'(r)EE J rfr'p^(r')pG(|r-r'|). (24) 

In almost all cases of interest we can picture V as aris- 
ing from a fixed, external charge distribution p^^^ , and we 
should apply the same separation of the Coulomb poten- 
tial to the fixed as well as mobile charges. This will allow 
us to express ([25]) in a suggestive and compact form using 
the sum of the fixed charge density and the equilibrated 
mobile charge density profile, Ptot(i") = P'^i^) + Plxti''^)- 
Thus we write 

V(r) = Vo(r) + Vi(r). (25) 

Just as vi (r) is the electrostatic potential due to a point 
charge convolved with a Gaussian of width cr, so Vi is the 
electrostatic potential due to the convolution of the fixed 
external charge density Pg^t ^it^ that same Gaussian as 
defined in ((TTI) and used in and ([M]). Given the 
definition of Vq and Vi , we write analogously 

VR(r) = Vo(r) + VRi(r), (26) 
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thus defining VRi(r), the slowly- varying component of 
the restructured electrostatic potential. 

Using (|23p . this slowly- varying component is exactly 
given by 



1 



(27) 



We see that Vri is the electrostatic potential arising 
from the Gaussian smoothing of the total charge density, 
P?ot(r). 

Equation (|77|) appears as the standard definition of the 
electrostatic potential due to charge density from first- 
year textbooks like [2^. Thus the LMF equation for 
the slowly- varying component of the restructured electro- 
static potential may be viewed as simply doing classical 
electrostatics on an equilibrium, smoothed total charge 
density and defining a single external restructured elec- 
trostatic potential. In fact we may in general represent 
the slowly- varying LMF equation (|27|) as a modified Pois- 
son's equation, 



V2VRi(r) 



47r 



PR,tot(r), 



(28) 



based on the smoothed total charge density. This point 
was first made for the uniformly-charged- wall model sys- 
tem [i^l , and it is immediately generalizable to the LMF 
equation dealt with here. The existence of this Poisson- 
like equation dictating a modified electrostatics empha- 
sizes that, in order for LMF theory to be valid for elec- 
trostatics, all charge-charge interactions, whether fixed, 
bound, or mobile, should be mapped via LMF theory 
to short-ranged interactions in the restructured electro- 
static potential Vr. In we emphasized the value 
of this modified-electrostatics perspective in physically 
interpreting molecular simulations and in understand- 
ing why and when spherical Gaussian truncations, using 
Vo{r) alone, fail in nonuniform environments. 



IV. FURTHER SIMPLIFICATIONS DUE TO 
CHARGE SMOOTHING 

Understanding the LMF equation as standard elec- 
trostatics for a Gaussian-smoothed charge density also 
allows for simpler solutions of the LMF equation than 
might be initially expected. As an example, in [Tl| . we 
were able to use a one-dimensional LMF equation to treat 
a system that has a charge density profile with an explic- 
itly three-dimensional character. In this section, we give 
the reasoning behind such simplifications and argue that 
they should be expected generally in many physically rel- 
evant cases. 



A. Simple Corrugated Surface 

In we treated SPC/E water confined to a slab 
by two empirical Pt(lll) surfaces [1^ that order the 
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FIG. 2: A (111) surface and the proximal charge density. 

(a) Sketch of the rectangular unit cell for a (111) surface show- 
ing the topmost three layers of atoms, labelled A, B, and C. 

(b) Corresponding p''{x,y) in eo/K^ for the 1.5 A layer closest 
to the wall. The density is projected onto the xy-plane, and 
the grayscale (color online) indicates the value of the charge 
at each point. Black dots indicate collected data points. 



water molecules at the surface, attracting the oxygen 
atoms to localized binding sites. Details of the simula- 
tion methods are available in [ll| . In brief, we simulated 
1054 SPC/E molecules confined between the two walls 
with 500 ps of equilibration and 1.5 ns of subsequent 
data collection. Simulations using slab-corrected three- 
dimensional Ewald sums. LMF theory with a — 4.5 A, 
and LMF theory with a = 6.0 A were carried out. Dur- 
ing those runs, the timestep chosen was 1.0 fs. However, 
for self-consistent solution of the LMF equation, we first 
simulated sets of 10 shorter runs with 25 ps of equili- 
bration and 50 ps of accumulation with a timestep of 
2.5 fs. These 10 shorter runs began with distinct initial 
conditions; we observed that averaging the p'^ resulting 
from such a set of simulations better isolated the system's 
equilibrium response to changes in Vr from the natural 
fluctuations in the charge density. For the larger a, con- 
vergence was achieved within 3 iterations, and for the 
smaller cr, 10 iterations were required for self-consistency. 
This approach certainly does not represent a numerically 
optimized solution of the LMF equation and is used here 
only to demonstrate the accuracy possible in principle 
when the full LMF theory is used. We will address the 
need for a more comprehensive and numerically efficient 
approach to the solution of the LMF equation and gen- 
eral timing issues further in section IVl 

In contrast to the smooth hydrophobic surface also 
dealt with in that paper, the empirical (111) surface was 
intended to represent detailed and specific ordering of 
water molecules at an atomically-corrugated surface. A 
schematic of the rectangular unit cell of a (111) surface 
is shown in figure 2(a) demonstrating the spacing of the 
three topmost layers of atoms. The empirical Pt(lll) 
potential [1^ that we use docs not have explicit atoms 
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or charges. Rather the potential simply encompasses the 
presence of binding sites through exponential functions 
and various cosine functions to represent the periodicity 
of the surface. 

As expected from the form of the potential, the charge 
density depends on x and y in addition to z. In 
figure 2(b) the charge density p'^{x,y) in the 1.5 A 
layer nearest the (lll)-wall is projected onto a unit cell. 
As previously demonstrated [2^, the surface induces a 
charge density profile very nonuniform in the x- and y- 
directions. The distinct periodicity is due to the fact that 
the spacing between binding sites (2.77 A) is nearly com- 
mensurate with the first peak in the bulk water gooi''') 
2.8 A). However despite a charge density that was ex- 
plicitly a function of r, we obtained accurate results for 
system properties using Vr{z) [ii|. Wc posited that this 
was the case due to the charge smoothing described in 
the previous section. 

Here we seek to demonstrate this more quantitatively 
by probing the charge density profile away from the sur- 
face within hexagonal prisms defined by three distinct 
hexagonal binding regions - A, B, and C - as labeled 
in figure |2(a) Figure [3]Ja) shows the p''{z) measured 



along each hexagonal prism and averaged laterally over 
the hexagonal area. These density profiles are quite dis- 
tinct both from each other and from the full laterally- 
smoothed p'^{z). 

The Gaussian smoothing of charge density inherent in 
the LMF treatment allows us to distinguish the coherent 
and incoherent variations in charge density seen in fig- 
ure [3](a). As shown in figure [3]Jb), convoluting /3^(r) with 
Gaussians in the x- and y-directions and subsequently 
analyzing sites A, B, and C leads to a density, Pjij' 
that is indistinguishable from the p'^{z) calculated ini- 
tially as only a function of z. This equivalence is due 
to that fact that the corrugations in effect induce local 
random noise in the density profile that is therefore av- 
eraged out by the Gaussians of width greater than the 
corrugation. In contrast, the presence of the confining 
surface leads to a p"^" density profile that docs not aver- 
age out in the z-direction as seen in figure[3J;. In essence, 
the LMF smoothing approach shows when variations in 
charge density profiles cancel, and when they do not. 



B. Simplifications Still Hold For Surfaces with 
Fixed Charges 
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FIG. 3: Effect of smoothing with a Gaussian of width a = 
4.5 A on p''(r) between two empirical (111) walls, (a) The 
charge density p'^i^z) for the three hexagonal prisms labeled 
A, B, and C in figure 2(a) (b) Gaussian smoothing of these 
density profiles in the x- and j/- direct ions, (c) Full three- 
dimensional Gaussian smoothing of the p'(z) profiles. The 
plots in (a) and (b) focus on data near the left wall, but the 
plot in (c) displays data for the full slab of water between the 
walls. 



charges associated with those molecules would also be 
included in the mobile p'. 

If, instead, the charges are held fixed, as is more likely 
for solid surfaces in classical simulations of solid-liquid in- 
terfaces, then (^ai^) = '/'ct,no(r) + 9aV(r). In such cases, 
we may no longer approximation VR,(r) as Vr(z), but 
a similar simplification for the smooth Vri will hold, 
as mi ght be expected. Such a potential for silica sur- 
faces [29| is currently being employed in our group to 
study the silica-acetonitrile interface [soj . 

For surfaces with fixed charges, the imposed external 
electrostatic potential V(r) would be 



However, a reasonable criticism of the previous surface 
is that it was not atomically realistic, even on a classi- 
cal level. More realistic metallic surfaces would require 
a treatment of image charges, which we will not address 
here. However, molecular surfaces modeled with stan- 
dard force fields would at least have charges at atomic 
sites near the surface. If the charges are free to move 
throughout the simulation as would generally be the case 
for biological membranes or liquid- liquid interfaces, ab- 
solutely no change in tactic would be necessary. The 



sites i 

where the prime indicates that these surface charges will 
need to be accounted for in the surface extending to in- 
finity in the x- and y-directions. The pair interactions 
representing the excluded volume of these atomic sites 
would be included in i^ine.Q- The external electrostatic 
potential V may be split into short-ranged component 
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FIG. 4: Vi{z) for charged (lll)-surface for two different a. 
This approximation to Vi (r) is nonneghgible within 2a of the 
surface. 

Vo and long-ranged component Vi as 

E 7^o(l'--^il)' (30) 

sites i 

Vi(r)= J2'^v,i\r~r,\). (31) 

sites i 

Now Vo is a simple minimum image sum over short- 
ranged pair interactions. The long-ranged component 
Vi still must include the effect of surface sites extending 
out to infinity in the x- and y-dircctions. but since Vi 
may be interpreted as the electrostatic potential due to a 
Gaussian smoothing of the fixed charge sites, we expect 
that 

Vi(r)~Vi(^). (32) 

The sole requirement would be that a should be larger 
than relevant lateral spacing of sites on the surface, such 
that the charge density will be sufficiently smoothed in 
the lateral directions. In such a case, Vr is a function of 
r, but all dependence on x and y is contained in Vo(r); 
and we would expect the self-consistent Vri(z) to hold 
to a very good approximation. 

As a quick demonstration of the validity of these 
ideas, we construct a (lll)-surfacc with lattice constant 
a = 3.92 A as for the model (lll)-surface discussed in 
section IIV Al However, we arbitrarily assign charges of 
-1-1 Co for the first layer of atoms, —2 eg for the second 
layer of atoms, and -|-1 Cq for the third layer of atoms. 

The corresponding Vi{z) would be determined by 
smoothing out the charges in each layer laterally, such 
that the charged (lll)-surface is now approximated by 
three uniformly charged planes. With this representa- 
tion, we may exactly express Vi{z) as the sum of three 
analytical functions, 

Vi(z)= ^G(z,z,), (33) 

?G layers 



where Xi is the surface charge density of each layer, and G 
is the one-dimensional Gaussian-smoothed Green's func- 
tion for a planar charge distribution [1, 0] analogous to 
vi{r) in Eq. (fT^ for a point charge. 

Vi{z) is plotted in figure|3]and is an important contrib- 
utor to the forces near the surface. However, since the 
surface is net neutral, it decays essentially to zero within 
a distance of approximately 2a from the surface. 

We may calculate both V(r) and Vo(r) easily by con- 
structing various configurations of one -1-1 unit test 
charge above this (lll)-surface in the dlpoly2.16 molec- 
ular dynamics package. We exactly calculate Vo(r) to 
within the simulation energy precision, by carrying out 
single point energy calculations for a +1 eo test charge 
above various (x, ?/)-sitcs interacting with each particle 
in the surface via vo{r). This particle is a test charge 
in the sense that it measures the Vq generated by the 
charges in the surface without disturbing the organiza- 
tion of charges within the surface. Thus the LMF esti- 
mation of the total V(r) is 

V(r) ~ VLMF(r) ~ Vo(r) + Vi(z), (34) 

with Vi(z) given analytically by ((33|) . 

In order to determine V(r) independently, we employ 
slab-corrected three-dimensional Ewald sums as a way 
to estimate that function (sij . Since Ewald sums re- 
quire net neutrality, we construct two surfaces with test 
charges ~ the surface we are interested in with a +1 
test charge above it and a second mirror surface with 
opposite charges with a —1 test charge above it. The 
two (lll)-surfaces have lateral size 30.49 A X 28.806 A 
with Lz ~ 280 A, and we use a = 0.34 A~^ and 
kmax = (13, 12, 120) for the corrected three-dimensional 
Ewald sums. We separate the two surfaces within the 
simulation cell as much as possible in order to decouple 
each surface and test-charge pair from the other pair. 

The energies yielded from these single-point-energy 
calculations J7Ewaid(r) allow us to estimate V(r), provided 
that we subtract the contribution from test charge-test 

charge interactions (z) and the constant contribution 

from intra-wall particle interactions (?7waii) as follows, 

V(r) ~ VEwald(r) = ^ {?/Ewald(r) - U+-{z) - C/walls} • 

(35) 

This is approximate since the Ewald sums will propagate 
the test charges laterally into all the periodic image cells, 
but given that the lateral area spans 11x6 repeats of the 
minimum rectangular (lll)-surface unit cell of 2.77 A x 
4.80 A, interactions between each surface and the lat- 
eral periodic images of the nearest test charge should be 
minimal. 

Shown in figure [S] is V{x,y,z) for z fixed at 2.26 A 
as treated by (a) Ewald sums, (b) the LMF approxima- 
tion with a = 4.5 A, and (c) the LMF approximation 
with a — 6.0 A. The similarity of the electrostatic po- 
tential surface for each a to the full Ewald calculation 
demonstrates that 4.5 A is greater than (Tmin for the fixed 
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FIG. 5: Electrostatic potential V(r) at a distance of 2.26 A 
above the charged surface as calculated by (a) corrected 
three-dimensional Ewald sums and by LMF estimation with 
(b) a = 4.5 A and (c) a = 6.0 A. LMF estimation entails 
approximating V(r) as Vo(r) + Vxiz). 



surface at least. This stands to reason since, within a 
plane, particles are within 2.77 A of each other and the 
planes of particles are separated by about 2.26 A. There- 
fore (T = 4.5 A is greater than relevant particle spac- 
ings. The agreement between VLMF(r) and VEwaid(r) is 
expected but numerically nontrivial since |Vi| ^ |V| at 
this distance above the surface. Therefore, achieving this 
degree of accuracy is a result of the accurate addition of 
two large-magnitude terms, Vo(r) and Vi(z), to yield a 
much smaller magnitude VlmfIi")- 

In order to demonstrate this agreement more broadly. 



we again study sites A, B, and C as shown in figure 2(a) 
Figure [6] shows the difference 



AV(r) = VE„aid(r) - VLMF(r), 



(36) 



laterally averaged above each hexagonal site. This vari- 
ation AV(r) is relatively small for a a of either 4.5 A or 
6.0 A since either a is greater than any relevant spacings 
in the surface. Compared to the huge discrepancies that 
would arise from using solely Vo(r) and neglecting the 
several Volt contribution of \>\[z), these differences are 
miniscule. At larger z, there are slightly larger variations, 
which we believe are numerical artifacts arising from our 
use of a sharp cutoff radius (10.5 A for cr = 4.5 A and 
13.5 A for a = 6.0 A) in Vo- Given the small scale on 
these graphs, even very tiny errors will visibly accumu- 
late. 

We fully expect that we should be able to treat Vr as 



VR(r)~Vo(r)+VRi(z). 



(37) 



In fact, we expect this approximation for Vri to be even 
better than for Vi since the mobile charge density is in- 
cluded and will act to neutralize the fixed charged sites to 
an extent. The validity of approximating the Gaussian- 
smoothed equilibrium charge density profile p'" as solely 
a function of z for non-electrostatic corrugated surface 
confinement discussed in section IIV Al should also hold 
for the mobile charge density profile in this model sys- 
tem with fixed point charges. Thus we conclude that 
this system should be treatable by LMF theory with the 
same one-dimensional equation as used in [111 ] with the 
sole difference that Vq is a function of r. 
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FIG. 6: Plots of AV(r) in Volt s aver aged above the sites 
A, B, and C as displayed in figure 2(a) for smoothing lengths 
a = 4.5 A and a = 6.0 A. Data for site B is displaced vertically 
by 0.04 Volts and data for site C by 0.08 Volts. In each case 
the value of Vi{z) is several volts, but the correction AV(r) 
is substantially smaller. 



V. CONCLUSIONS AND OUTLOOK 

Based on the success of LMF theory in treatinga vari- 
ety of charged and molecular systems 0, H, H, [lOl llll. , 
we believe that LMF theory is a promising new approach 
that permits simple minimum-image simulations of 1/r 
interactions while still accurately accounting for the net 
additive long-ranged forces using the restructured exter- 
nal electrostatic potential Vr. While the systems treated 
thus far have been simple ionic and molecular systems, 
the LMF framework is much more broadly applicable 
to simulation systems of biomolecular and experimental 
interest. Current research using the LMF approach in 
the Weeks research group includes the collapse of model 
charged polymers 13211 and the behavior of acetonitrile 
near silica surfaces [30[ . 

Here we have sought to present the statistical mechan- 
ical foundations of the LMF treatment of charged sys- 
tems. We have highlighted the necessary but physically 
reasonable approximations leading to the LMF equation, 
and we have also derived the simple LMF equation for the 
electrostatic potential valid when treating point-charge 
models, both ionic and molecular. 

Furthermore, we have examined the useful conse- 
quences of understanding the smoothly-varying LMF 
electrostatic potential Vri as resulting from a Gaussian- 
smoothed charge density. In a previous paper [ll[, we 
showed that this charge density can be physically enlight- 
ening. Here we instead show the practical utility of this 
smoothing. The convolution of the charge density with 



11 



a Gaussian of width a allows us to solve much simpler 
LMF equations than might originally be believed. For 
confinement into a slab geometry, regardless of whether 
the confinement has nonelectrostatic or electrostatic cor- 
rugations, the relevant LMF equation will depend only 
on z to a very good approximation. 

In future work, wc will examine the application of LMF 
theory to bulk site-site fluids in much more detail [2^ and 
also demonstrate very simple corrections that allow much 
more accurate thermodynamics from GT models j33| . 
Another area which we hope to explore in greater detail is 
methods for optimally solving the LMF equation during 
simulation. As alluded to in [ll| , the writing of the LMF 
equation as a Poisson equation with a Gaussian smoothed 
charge density suggests that connection with work on ef- 
ficient Poisson solvers would be numerically helpful, and 
that siinple approximations may prove surprisingly ac- 
curate [lO|. Furthermore, since Poisson's equation may 
be written as the minimum of a functional [3J|, using a 
modified Car-Parinello scheme [s^ to slowly evolve to 
the averaged, equilibrium solution of the LMF equation 
seems like another fruitful path to pursue. In our (un- 
optimizcd) simulations of LMF systems, we observed a 
speed-up of at least a factor of four relative to Ewald 
summation. However, such timing questions need to be 
studied in much greater detail once an optimized path to 
LMF solution is developed. 

Regardless of such possibilities for solution of the LMF 
equation, LMF theory is still an inherently equilibrium 
theory. Given the success of a scaling principle of dynam- 
ics that states that equilibrium, dynamic properties will 
be well-captured in a reference system that reproduced 
static, structural properties [s^, many dynamical prop- 
erties for equilibrium systems could possibly be modeled 
quite well in our LMF-mapped systems. But such a con- 
jecture needs to be tested, and, regardless, much further 
theoretical work is required to develop the analog of LMF 
theory for dynamically-evolving systems. 
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APPENDIX: RELATION BETWEEN THE GT 
voir) AND REACTION FIELD TRUNCATIONS 



The standard short-ranged RF pair interaction 
t'o,RF(|rn — Ttl) is proportional to the potential energy be- 
tween a unit point charge at and a unit point charge at 
r„ surrounded by a uniform neutralizing spherical charge 
distribution terminating at the sharp cutoff radius Rc 
of the RF model. Anticipating that we may want a 
more symmetric description of the charge distributions 
at Tt and r„, we can also relate the RF wq.rf to the dif- 
ference between the potential energy of two unit point 
charges at and r„ and the potential energy of a unit 
uniform spherical charge distribution at r„ and a unit 
point charge at rt. But the subtracted term is not sym- 
metric and the smoother charged-clouds interaction [l3| 
''^0,cc(^) can be found by subtracting the potential en- 
ergy of two unit uniform spherical charge distributions 
at r„ and rj (with cutoff radius i?c/2) from the potential 
energy of two unit point charges at and r^. The func- 
tional form wo,cc('') was previously arrived at via a less 
symmetric argument [14| . The functional form VQ,cc{f) 
is different than that of wo,rf('')- 



In contrast, the z;o(|r„-rt|) = erfc(|r„-rt|/(T)/|r„-rt| 
chosen for LMF theory can be related to either the energy 
between a unit point charge at and a unit point charge 
at Tn surrounded by a neutralizing Gaussian charge dis- 
tribution of width a or equivalently to the difference be- 
tween the potential energy of two unit point charges at 
rt and r„ and the potential energy of two unit Gaus- 
sian charge distributions of width a j \f2 at r,i and rt 0]. 
The same functional form v^ir) = crk{r/a)/r results in 
either case. 
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